#include <gsl/gsl_rng.h>
#include <math.h>
#include <stdio.h>

int Nmesh;
unsigned int *seedtable;
gsl_rng *random_generator;
double amp_out[10240],phase_out[10240];

#define  PI          3.14159265358979323846
#define  ROOT2       1.414213562

void get_ij_k_rands23_(long long int *, long long int *, long long int *, 
		     float *, float *);

void initialise_random_seed_(int *, long long int *);
void initialise_random_seed_mill_(int *, long long int *);
void set_standard_table(unsigned int *, gsl_rng * );
void set_millennium_table(unsigned int *, gsl_rng *);
void principal_(int *, int *, int *, int *);
int FatalError(int );



void get_ij_k_rands23_(long long int *nft, long long int *ii, long long int *jj, 
		     float *ggvar1, float *ggvar2)
{
  int k,veclen;
  int ix0,iy0;
  int ix,iy,iz;
  int ixx,iyy;
  veclen = *nft/2;
  double ampl,phase;
  float multi;
  float g1,g2;
  int isign;
  int iflag;

  ix = *ii;
  iy = *jj;
  ix0 = ix;
  iy0 = iy;
  ixx = (*nft-ix)%(*nft);
  iyy = (*nft-iy)%(*nft);

  gsl_rng_set(random_generator, seedtable[ix * (*nft) + iy]);

  if (ix>veclen) ix=ix-*nft;
  if (iy>veclen) iy=iy-*nft;



  for (k=0; k<veclen;k++){
    iz = k;
    if (iz>veclen) iz=iz-*nft;
    phase = gsl_rng_uniform(random_generator) * 2 * PI;
    do
     ampl = gsl_rng_uniform(random_generator);
     while(ampl == 0); 
 
    amp_out[k] = ampl;
    phase_out[k] = phase;
       
  }


  /* May Now generate complementary k=0 value */

iflag = 0;   

     /* Now decide which k=0 mode to use */

 if (ix0==0){
       if (iy0>=Nmesh/2) {
	 iflag = 1;
      }} else {
	   if (ix0>=Nmesh/2){
	     iflag = 1;
	   }
 }

  /* iflag==1 requires generation of second ampl/phase pair */


if (iflag == 1){
         gsl_rng_set(random_generator, seedtable[ixx * (*nft) + iyy]);
         phase = gsl_rng_uniform(random_generator) * 2 * PI;
         do
           ampl = gsl_rng_uniform(random_generator);
        while(ampl == 0);
        amp_out[0]   =  ampl;
        phase_out[0] = -phase;
   }

 for (k=0; k<veclen;k++){
    principal_(&ix,&iy,&iz,&isign);
    phase = phase_out[k];
    ampl  = amp_out[k];


     phase = -phase;  

     if ((ix>0)&&(iz<0)) phase = -phase;
     if ((ix<0)&&(iz>0)) phase = -phase;

     if ((ix==0)&&(  ((iy>0)&&(iz<0)) || ((iy<0)&&(iz>0)))) phase = -phase;


     if (ampl!=0){
       g1 = sqrt(-log(ampl)) *ROOT2 * cos(phase);
       g2 = sqrt(-log(ampl)) *ROOT2 * sin(phase) * (float)isign; 
     }else{
       g1 = 0;
       g2 = 0;
     }


     ggvar1[k] = g1;
     ggvar2[k] = g2;

 }
  	 
}
















void initialise_random_seed_(int *Seed, long long int *nft){
  Nmesh = *nft;
random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
gsl_rng_set(random_generator, *Seed);
  if(!(seedtable = malloc(Nmesh * Nmesh * sizeof(unsigned int))))
    FatalError(4);
  set_standard_table(seedtable,random_generator);
}

void initialise_random_seed_mill_(int *Seed, long long int *nft){
  Nmesh = *nft;
random_generator = gsl_rng_alloc(gsl_rng_ranlxd1);
gsl_rng_set(random_generator, *Seed);
  if(!(seedtable = malloc(Nmesh * Nmesh * sizeof(unsigned int))))
    FatalError(4);
  set_millennium_table(seedtable,random_generator);
}

void set_standard_table(unsigned int *seedtable, gsl_rng * random_generator)
{
  int i,j;
                        /*  Standard VRS WhichRandMethod == 1 */
  for(i = 0; i < Nmesh / 2; i++)
    {
      for(j = 0; j < i; j++)
	seedtable[i * Nmesh + j] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedtable[j * Nmesh + i] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i; j++)
	seedtable[(Nmesh - 1 - i) * Nmesh + j] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedtable[(Nmesh - 1 - j) * Nmesh + i] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i; j++)
	seedtable[i * Nmesh + (Nmesh - 1 - j)] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedtable[j * Nmesh + (Nmesh - 1 - i)] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i; j++)
	seedtable[(Nmesh - 1 - i) * Nmesh + (Nmesh - 1 - j)] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedtable[(Nmesh - 1 - j) * Nmesh + (Nmesh - 1 - i)] = 0x7fffffff * gsl_rng_uniform(random_generator);
    }

}






void set_millennium_table(unsigned int *seedtable, gsl_rng * random_generator)
{
  unsigned int *seedmill;
  int Nmill, i, j;

  Nmill = 2560;

  if(!(seedmill = malloc(Nmill * Nmill * sizeof(unsigned int))))
    FatalError(4);

  if(Nmesh > Nmill)
    FatalError(5);

  for(i = 0; i < Nmill; i++)
    {
      for(j = 0; j < i; j++)
	seedmill[i * Nmill + j] = 0x7fffffff * gsl_rng_uniform(random_generator);

      for(j = 0; j < i + 1; j++)
	seedmill[j * Nmill + i] = 0x7fffffff * gsl_rng_uniform(random_generator);
    }


  for(i = 0; i < Nmesh / 2; i++)
    {
      for(j = 0; j < i; j++)
	seedtable[i * Nmesh + j] = seedmill[i * Nmill + j];

      for(j = 0; j < i + 1; j++)
	seedtable[j * Nmesh + i] = seedmill[j * Nmill + i];

      for(j = 0; j < i; j++)
	seedtable[(Nmesh - 1 - i) * Nmesh + j] = seedmill[(Nmill - 1 - i) * Nmill + j];

      for(j = 0; j < i + 1; j++)
	seedtable[(Nmesh - 1 - j) * Nmesh + i] = seedmill[(Nmill - 1 - j) * Nmill + i];

      for(j = 0; j < i; j++)
	seedtable[i * Nmesh + (Nmesh - 1 - j)] = seedmill[i * Nmill + (Nmill - 1 - j)];

      for(j = 0; j < i + 1; j++)
	seedtable[j * Nmesh + (Nmesh - 1 - i)] = seedmill[j * Nmill + (Nmill - 1 - i)];

      for(j = 0; j < i; j++)
	seedtable[(Nmesh - 1 - i) * Nmesh + (Nmesh - 1 - j)] =
	  seedmill[(Nmill - 1 - i) * Nmill + (Nmill - 1 - j)];

      for(j = 0; j < i + 1; j++)
	seedtable[(Nmesh - 1 - j) * Nmesh + (Nmesh - 1 - i)] =
	  seedmill[(Nmill - 1 - j) * Nmill + (Nmill - 1 - i)];
    }

  free(seedmill);
}

int FatalError(int errnum)
{
  printf("FatalError called with number=%d\n", errnum);
  fflush(stdout);
  abort();
}

